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Theory of dilution effect in orbital ordered system is presented. The e g orbital model without spin degree 
of freedom and the spin-orbital coupled model in a three-dimensional simple-cubic lattice are analyzed by the 
Monte-Carlo simulation and the cluster expansion method. In the e g orbital model without spin degree of free- 
dom, reduction of the orbital ordering temperature due to dilution is steeper than that in the dilute magnet. This 
is attributed to a modification of the orbital wave-function around vacant sites. In the spin-orbital coupled model, 
it is found that magnetic structure is changed from the A-type antiferromagnetic order into the ferromagnetic 
one. Orbital dependent exchange interaction and a sign change of this interaction around vacant sites bring 
about this novel phenomena. Present results explain the recent experiments in transition-metal compounds with 
orbital dilution. 
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I. INTRODUCTION 

Impurity effect in correlated electron system is one of the 
attractive themes in recent solid state physicsJ* 2 - The well 
known example is doping of non-magnetic impurity in high 
Tc superconducting cuprates; a small amount of substitu- 
tion of Cu by Zn dramatically destroys the superconductivity. 
Non-magnetic impurity effect in the low-dimensional gapped 
spin system is another example. A few percent doping of Zn 
or Mg, which does not have a magnetic moment, into two- 
leg ladder systems, e.g. SrCu2C»3, and spin-Peierls systems, 
e.g. CuGeC»3, induces long-range orders of antiferromag- 
netism (AFM) i 3 i 4 i 5 ' 6 ' 7 Impurity effect in charge and orbital 
ordered state is also studied in the colossal magnetoresistive 
manganites^ It is reported in a so-called half-doped mangan- 
ite Lao.5Cao.5Mn03 that a few percent substitution of Mn by 
Cr collapses the charge/orbital order associated with the AFM 
one and induces a ferromagnetic metallic state. Because of no 
e g electrons in Cr 3+ , unlike Mn 3+ with one e g electron, Cr is 
regarded as an impurity without orbital degree of freedom. 

Recently, impurity doping effect in an orbital ordered state 
is examined experimentally in a more ideal material. Mu- 
rakami et al. have studied substitution effect in an orbital 
ordered Mott insulator KCUF3 with the three dimensional 
(3D) Perovskite crystal structured A Cu 2+ ion in the cubic- 
crystalline field shows the (t2 g ) 6 (e g ) 3 electron configuration 
where one hole has the orbital degree of freedom. The long- 
range orbital order (OO), where the d i_ z i- and c/j_ v 2-like 
orbitals are aligned with a momentum (n, 7t, n), was observed 
at room temperatures by several experiments. Since the AFM 
spin ordering temperature is 39K, which is much lower than 
the OO temperature (> 1200K), a substitution of Cu by Zn, 
which has an electron configuration (t2g) 6 (eg) 4 , is regarded as 
an orbital dilution. It was revealed by the resonant x-ray scat- 
tering experiments in KCui_ r Zn v F3 that the OO temperature 
decreases with doping of Zn monotonically and the diffraction 
intensity at (3/2 3/2 3/2) disappears around x = 0.45. At the 
same Zn concentration, the crystal symmetry is changed from 
the tetragonal to the cubic one. That is to say, the OO disap- 
pears around x = 0.45. In dilute magnets, e.g. KMni_ r Mg v F3, 
the x dependence of the magnetic ordering temperature as well 



as the critical concentration where the magnetic order van- 
ishes are well explained by the percolation theory. 11 ' 12 On the 
contrary, the critical concentration in KCui_ v Zn v F3, where 
the OO disappears, is much smaller than the site-percolation 
threshold in a 3D simple cubic lattice, x p = 0.69. These ex- 
perimental observations imply that the dilute OO may belong 
to a new class of diluted systems beyond the conventional per- 
colation theory. 

Dilution effect in orbital ordered state was also examined 
experimentally in a mother compound of the colossal magne- 
toresisitive manganites, LaMn03. The long-range OO, where 
the dj x 2_ r 2- and c/ 3v 2_ r 2-like orbitals align with a momen- 
tum (n,n,0), appears below 780K. The A-type AFM order, 
where spins are aligned ferromagnetically in the xy plane 
and are antiferromagnetically along the z axis, is realized 
at 140K. Substitution of Mn 3+ by Ga 3+ , which has a 3<f 10 
electron configuration, corresponds to both the orbital and 
spin dilution . 13 ' 14 i 15 ' 16 ' 17 i 18 From the x-ray diffraction and X- 
ray absorption near-edge structure (XANES) experiments, the 
tetragonally distorted MnOg octahedra become regular cubic 
ones around the Ga concentration x = 0.6. That is, the OO 
disappear around x = 0.6 which is smaller than the percola- 
tion threshold x p = 0.69 for the simple cubic lattice. Differ- 
ence between LaMni_ x Ga x 03 and KCui_ x Zn x F3 is seen in 
the magnetic structure. Blasco et al. observed by the neu- 
tron diffraction experiments in LaMni_ v Ga r 03 that the fer- 
romagnetic (FM) component appears by substitution by Ga 
and increases up to x = 0.5. This change of the magnetic 
structure from the A-type AFM to FM was also confirmed by 
the magnetization measurements. This FM component can- 
not be attributed to the itinerant electrons through the double 
exchange interaction, since the electrical resistivity increases 
with increasing x. These phenomena are in contrast to the 
conventional dilute magnets where the ordering temperature 
is reduced, but the magnetic structure is not changed. Far- 
rell and Gehring presented a phenomenological theory for the 
magnetism in LaMni_j-Ga r 03 M They noticed that a volume 
in a GaOe octahedron is smaller than that in a MnOe- Un- 
der an assumption that the Mn 3d orbitals around a doped Ga 
tend to be toward the Ga, the magnetic structure change was 
examined. 
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In this paper, a microscopic theory of dilution effects in the 
e g orbital degenerate system is presented. We study the dilu- 
tion effects in the e ? -orbital Hamiltonian without the spin de- 
gree of freedom, termed Mj [see Eq. (TTOlil. and the spin and 
e g orbital coupled one, termed J$?st [see Eq. (|9)]. The classi- 
cal Monte-Carlo (MC) method in a finite size cluster, as well 
as the cluster expansion (CE) method is utilized. It is known 
that, in the classical ground state of J%t without impurity, a 
macroscopic number of orbital states are degenerated due to 
frustrated nature of the orbital interaction. We demonstrate 
numerically that this degeneracy is lifted at finite temperature. 
It is shown that the OO temperature decreases rapidly with 
increasing dilution. From the system size dependence of the 
orbital correlation function in the MC method, the OO is not 
realized at the impurity concentration x = 0.2. The results ob- 
tained by the CE method also show rapid quenching of OO by 
dilution in comparison with dilute spin models. These results 
are interpreted that orbitals around impurity sites are changed 
so as to gain the remaining bond energy. This is a conse- 
quence of the bond-direction dependent interaction between 
the inter-site orbitals. In the analyses of the spin-orbital cou- 
pled model, it is shown that the A-type AFM structure realized 
in x = is changed into FM one by dilution. This is explained 
by changing a sign of the magnetic exchange interaction due 
to the orbital modification around impurity sites. Implications 
of the present microscopic theory and the experimental results 
in KCui_ A Zn x F3 and LaMni_ x Ga T 03 are discussed. 

In Sect. II, the model Hamiltonian for the e g orbital degree 
of freedom in a cubic lattice and the spin-orbital coupled one 
are introduced. In Sect. Ill, the classical MC simulation and 
the CE method are presented. Results of the numerical analy- 
ses in 3%j and ,J^st are presented in Sects. IV and V, respec- 
tively. Section VI is devoted to summary and discussion. A 
part of the numerical results for the e g orbital model have been 
briefly presented in Ref.[l9l 



II. MODEL 

Doubly degenerate e g orbital degree of freedom is treated 
by the pseudo-spin (PS) operator with magnitude of 1/2. This 
operator is defined by 

T i = \ £ d t rs t V d iy's, (1) 

syf 

where djy s is the annihilation operator of an electron with spin 
s(=T, I) and orbital y(= 3z 2 — r 2 7 x 2 — y 2 ) at site i, and a are 
the Pauli matrices. Occupied orbital is represented by an angle 
6 of PS. The eigen state of the z-component of PS with an 
angle 8 is 




For example, = 0, 2n/3, and 4n/3 correspond to the states 
where the d 3z 2_ r 2, d 3y 2_ r 2, and d 3x 2_ r 2 orbitals are occupied 
by an electron, respectively. It is convenient to introduce the 



linear combinations of the PS operators defined by 

Tf=cosl — J5f-sinl — \Tf, (3) 

with / = (x, y, z) and a numerical factor (n x ,n y ,n z ) = (1,2,3). 
These are the eigen operators for the dy2_ r 2 orbitals. 

It is known that dominant orbital interactions in transition- 
metal compounds are the electronic exchange interaction and 
phononic one. The former is derived from the generalized 
Hubbard-type model with the doubly degenerate e g orbitals; 

^ele = £ (tj/dJ ys d J y s + H.c}j+UY j n, rl n, n 

(ij)yfs ij 

+ \u' £ n ir n iY + -K £ dj rs dl ys ,d irs ,d lY J4) 

i y^Y i y^Y ss ' 

where n,y = Y<s n iys — H s d}y s diy s - We define the electron trans- 
fer integral tjf between the a pair of the nearest neighbor- 
ing (NN) sites. The intra-orbital Coulomb interaction U, the 
inter-orbital one U', and the Hund coupling K. Through the 
perturbational expansion with respect to the NN transfer in- 
tegral under the strong Coulomb interaction, the spin-orbital 
superexchange model is obtained. 20 ' 21 By assuming a relation 
U = U' + K, for simplicity, it is given as 




-^l(j-SrS j )(^ + T/tj + T/ + Tj), (5) 

where Sj is the spin operator at site i with a mgnitude of 
1/2, and I represents a bond direction connecting sites i and 
j. Amplitudes of the superexchange interactions are given as 
Ji [= t 2 /(U - 3K)] and J 2 {= t 2 /U) where t is the transfer in- 
tegral between the d^2_ r 2 orbitals along the z direction. 

The phononic interaction between the orbitals is derived 
from the orbital-lattice coupled model given by 

^rr = -gjrZQTTr 

im 

where a subscript m takes x and z. The first term represents the 
Jahn-Teller (JT) coupling with a coupling constant gjj. Two 
distortion modes in a 06 octahedron with the E g symmetry is 
denoted by Qf and QJ. The second term is for the JT phonon 
where and are the phonon coordinate and momentum, 
respectively, and is the phonon frequency. Subscripts k 
and 4 are the momentum and the phonon mode, respectively. 
Here, the spring constant between the NN metal and oxygen 
ions are taken into account. The interaction between orbitals 
and the uniform strain and the strain-energy, which are neces- 
sary in study of the cooperative JT effect, are not shown, for 
simplicity, in this equation. For convenience, the first and sec- 
ond terms in Eq.|6]are denoted by .^rb-iatt an d ^fatt> respec- 
tively. By introducing the canonical transformation defined 
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by 



<M, - c iki; 



(7) 



and neglecting the non-commutability between J#\ at and , 
the orbital and lattice degrees of freedom are decoupled 

„„22.23.24,25 



■sty 



JT 



(ij) 



latt- 



(8) 



The first term in this equation gives the inter-site orbital in- 
teraction with a coupling constant g = gj T /(3Ks) where Kg 
is a spring constant, and M\ att is given by the second term in 
Eq. (O, i.e. J%\att, where the phonon coordinate and momen- 
tum are replaced by and its canonical conjugate momen- 
tum p k e , respectively. 

The model Hamiltonian studied in the present paper is 
given by a sum of the above two contributions. Quenched 
impurity without spin and orbital degrees of freedom is de- 
noted by a parameter e, which takes zero (one), when site i 
is occupied (unoccupied) by an impurity. The Hamiltonian is 



given as 



.jfr ST = — 



(ij) V4 
(ij) V 



A 1 J 



- + T-t' + T- + T- 

A ' J ' J 



[ij) 



(9) 



Numerical results in this Hamiltonian is presented in Sect. [V] 
We also study dilution effect in the orbital model without spin 
degree of freedom. This model is given by taking Sj • Sj in 
Eq. (0 to be zero. This procedure may be justified in the 
diluted orbital system of KCui_ v Zn v F3 where the Neel tem- 
perature (Tn) is much below the OO temperature Too The 
explicit form of the e g orbital model without spin degree of 
freedom is given by 



■SIT 



(ij) 



(10) 



where J(— 2g + 3/i/4 — 72/4) is the effective coupling con- 
stant. Numerical results of this model Hamiltonian are pre- 
sented in Sect.HVl 



III. METHOD 

In order to analyze the model Hamiltonian introduced 
above by using the unbiased method, we adopt mainly the 
classical MC simulation in finite size clusters. The orbital PS 
operator is treated as a classical vector defined in the T z — T x 
plane, i.e. T- z = (l/2)cos0; and Tf = (l/2)sin0; where 0, is 
a continuous variable. As well as the conventional Metropolis 
algorithm, the Wang-Landau (WL) method is utilized^ This 



is suitable for the present spin-orbital coupled model where 
the energy scales of the two degrees are much different with 
each other. In order to calculate the density of state, g(E), with 
high accuracy in the WL method, we take that the minimum 
energy edge E m i n in g(E) is higher a little than the ground 
state energy Eqs, an d assume g(£bs < E < E m i„) =0. As a 
result, the present MC simulation is valid above a character- 
istic temperature r m ; n which is determined by \E m i„ — Eqs I ■ 
This situation will be discussed in Sect. |IV] in more detail. 
The simulations have been performed in L x L x L cubic lat- 
tices (L = 12 ~ 18) with the periodic-boundary condition. In 
the Metropolis method, for each sample, 3 x 10 4 — 1 x 10 5 MC 
steps are spent for measurement after 8x 10 3 — 2 x 10 4 MC 
steps for thermalization. Physical quantities are averaged over 
20 — 80 samples at each parameter set. In the WL method, the 
final modification factor— is set to be //;„„/ = exp(2~ 27 ). Af- 
ter calculating the density of states, 2x 10 7 MC steps are spent 
for measurement. 

To supplement the classical MC simulation, the ordering 
temperatures are also calculated by utilizing the CE method. 
We apply the CE method proposed in Ref. |27| to the present 
orbital model. For a given impurities configuration {e} in a 
lattice with sites, the OO parameter is given as 



M {e} =^^^1 PH 



with the density matrix 

Pn{e} 



Tr N e~ 



(11) 



(12) 



where Tr^ represents the trace over the PS operator at sites 
with £,■ = 1 in a crystal lattice, and J^ £ j is the Hamiltonian 
with impurity configuration {e}. The OO parameter per site 
is obtained by averaging about all possible impurity configu- 
ration {e} as 



M = 



1 



(l-x)N 



< M W> {e} - 



(13) 



where x is the impurity concentration. In the CE method, a 
cluster consisting of m sites, termed {m}, is considered, and 
the PS operators which do not belong to {m} are replaced by 
stochastic variables <7,-. Here we take (T : X ,T?) = (0,(7,). The 
effective Hamiltonian thus obtained is denoted as J#{ e y{ m }{ a } 
where {a} is a set of <7,, and the corresponding density matrix 
is 



P{e}{m}{a} 



exp(-/3^ {£}{m}{cr} ) 
Tr {m} exp (-PJf{e}{m}{o}) 



(14) 



where Tr{ m j represents the trace over the PS operators in a 
cluster {m}. We expand M{ £ j into a series of cluster averages 
as follows, 



^}=IIII(-1)*-" 

m=\ {m}k=l{k} 



x Tr 




£p{e}{*}{tT} 



(15) 





X P 



FIG. 1: (a) Scematic picture of the degenerate PS configurations 
termed the type-(I) degeneracy, and (b) that of the type-(II) one. 



where E{m} i s taken over all possible clusters consisting of m 
sites, and is taken over all subclusters of k sites belong- 
ing to a given {m}. The variable <7, takes 1/2 or —1/2 by a 
probability of 



1+2M 



1-2M 



(16) 



By solving Eqs. (fT3])-(fT6l) self-consistently, the order param- 
eter and the ordering temperature are obtained as a function 
of impurity concentration. In the present study, we adopt the 
CE method in the two-site cluster approximation, i.e. m = 2. 
It was shown that, even in the two-site cluster approximation, 
the obtained results show good accuracy in the case of the fer- 
romagnetic Heisenberg model in a simple cubic lattice; devia- 
tions from the results by other reliable methods are about 2% 
for the critical impurity concentration!^ To compare the re- 
sults in the classical MC simulation, the ordering temperature 
is also calculated in the classical version of the CE method 
where the traces in Eqs. ( fT4b and (|T~5T > are replaced by inte- 
grals with respect to the continuous variable Tf between 1 /2 
and -1/2. 



IV. DILUTION IN THE e g ORBITAL MODEL 

In this section, numerical results for dilution effects in the 
e g orbital Hamiltonian ( fTOb are presented. First, we show 
the results in the MC simulation without impurity. It is 
known that there is a macroscopic degeneracy in the mean- 
field (MF) ground state in Mr without impurity. 2 ^ This degen- 
eracy is classified into the following two types: (I) Consider 
a staggered-type OO with two sublattices, termed A and B, 
and momentum Q = (n, n, n). In the MF ground state, the PS 
angles in the sublattices are given by (GU, 9b)=(9, 9 + n) with 
any value of 9. Such continuous rotational symmetry is unex- 
pected from the Hamiltonian Mr where any continuous sym- 
metries do not exist. (II) Consider an OO with Q = {%, n, n) 
and (0/i, 9b)=(9o, 9q + n), and focus on one direction in three- 
dimensional simple-cubic lattice, e.g., the z direction. The 
MF energy is preserved by changing all PS in each layer 
perpendicular to the z axis independently as (9q, 9q + n) — > 




FIG. 2: (a) System size dependence of the orbital correlation function 
Moo { x = 0), and (b) that of the orbital angle function M am (x = 0). 
The minimum energy E m \ a in the WL method is taken to be 0.95£gs 
in (a) and 0.98£ G S in (b). 



(—9o, — 9q — n). These are schematically shown in Fig. Q] 
Both types of degeneracy are understood from the momentum 
representation of the orbital interaction, 



^=2/2>k£(k)Vk, 

k 



(17) 



with y/k = Pk i T£] and the 2x2 matrix E(k). By diagonaliz- 
ing £(k), we obtain the eigen values 



£±(k) 



± 



-ch 



CxCy CyC z 



(18) 



where c/ = cos ale/ with a lattice constant a. The lower eigen 
value /-(k) has its minima along (n,7Z,n) — (0,7T,7r) and 
other two-equivalent directions.— At the point F, the two 
eigen values E + (k) and £_ (k) are degenerate. That is, the or- 
bital states corresponding to these momenta are energetically 
degenerate in the MF level. A lifting of this degeneracy in the 
MF ground state has been examined from the view points of 
the order-by-disorder mechanism by utilizing the spin wave 
analyses i 28 i 30 ' 31 

Here we demonstrate the degeneracy lifting and appearance 
of the long-rage OO by the MC method. We introduce, for im- 
purity concentration x, the staggered orbital correlation func- 
tion 



^ooW = ^/{i(-iy e! T i 



1/2 



(19) 
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and the angle correlation function 

Mang(x) 



0.35 



2\ 1/2 



(20) 



where 



represents the MC average and N = L . The or- 



bital correlation at the momentum Q = (it, 7T, it) is represented 
by Moo(i), and the angle correlation M ang (x) takes one, when 
the orbital PS angle is 2itn /3 with an integer number n. There- 
fore, Mqq(x) and M ang (jc) are utilized as monitors for lifting 
of the type-(II) and (I) degeneracies, respectively. Temper- 
ature dependences of Mqq(x = 0) for various L are shown 
in Fig. 12 a). With decreasing temperature, calculated results 
for all L show a sharp increasing around T /J = 0.35. This 
increasing becomes sharper with the system size L. Below 
T/J = 0.08, Mqq(x = 0) takes a temperature-independent 
value of about 0.47. This flat behavior is attributed to the 
lowest energy edge E m \ n for the density of state calculated in 
the WL method, as explained in Sect. [HI] An extrapolated 
value of Mqo(x = 0) toward T — is close to 0.5 which in- 
dicates that the type-(II) degeneracy is lifted and the OO with 
the momentum Q = (it, it, it) is realized. Temperature depen- 
dences of M ang (x = 0) presented in Fig. |2jb) increase mono- 
tonically toward one in the low temperature limit. Almost no- 
size dependence is seen in M wg (x = 0). Therefore, the type- 
(I) degeneracy is also lifted and the PS angle is fixed. Both 
results indicate the long-range OO where the momentum is 
Q = (it, it, it), and the PS angles are (6a, 9b) = (9o, Oo + it) 
with O = 2nn/3. 

The temperature at which Moo(x = 0) and M ang (x = 0) 
change abruptly is around T/J = 0.33 corresponding to the 
OO temperature Too(x = 0). In more detail, this tempera- 
ture is determined by the finite-size scaling for the correlation 
length. This is calculated by the second-moment method; 



«to = 



1 



2sin(afc min /2) 



lMoo(x) 2 -M krm (x)^ (21) 



with 



M k . (x) 



M k . (x) 2 

"mm V / 



^ e /(Q-k)r i£ . T .| 



2v 1/2 



, (22) 



N(l -x) 

where kmm = (2it /L,Q,Q). The scaling relation for <fj (x) is 

$(x)=Lf[l 1 / v {T-T o(x)}}, (23) 

where v is the critical exponent for correlation length, and F 
is the scaling function. The correlation lengths § (x = 0)/L for 
various sizes cross with each other at Too(x = 0). In Fig. [3] 
we plot B,(x = Q)/L as a function of L}I V [T - T QO (x = 0)]. 
The scaling analyses work quite well for L = 10, 12, and 14. 
The OO temperature Too(x = 0) and the critical exponent v 
are determined by the least-square fitting for the polynomial 
expansion. We obtain as T 00 (x = 0)/J = 0.344 ±0.002 and 
v = 0.69 — 0.81, although statistical errors are not enough to 
obtain the precise value of v. 




L 1/v (T-T 00 )/Too 



FIG. 3: Scaling plot of the correlation length %(x = 0) for the 
staggered orbital correlation. Numerical data are obtained by the 
Metropolis algorithm. 
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FIG. 4: Impurity concentration dependence of the staggered orbital 
correlation function Mqq(x). System size is taken to be L = 18. Nu- 
merical data are obtained by the Metropolis algorithm. 



Now, we examine impurity effect in the OO. In Fig. [4] 
we present the staggered orbital correlation function Moo (x) 
for several impurity concentration x. Numerical data are ob- 
tained by the Metropolis algorithm in the classical MC method 
and the system size is chosen to be L = 18. First, we fo- 
cus on the region of x < 0.15. As shown above, Moo( x = 0) 
abruptly increases at Too(x = 0) ~ 0.347 and is saturated to 
0.5 in the low temperature limit. By introducing impurity, 
Moo(x > 0) does not reach 0.5 even at T /J = 0.01, and its 
saturated value in low temperatures gradually decreases with 
increasing x. Although the system sizes are not sufficient to 
estimate Mqo(x) in the thermodynamic limit, Mqq(x > 0) at 
zero temperature does not show the smooth convergence to 
0.5 in contrast to the diluted spin models. Beyond x = 0.15, 
results are different qualitatively; although Mqq(x) starts to 
increase around a certain temperature (e.g. T/J ~ 0.24 at 
x = 0.2), saturated values of Moo(x) in the low temperature 
limit are rather small. In order to compare the size depen- 
dences of Moo (x) < temperature dependences of Moo (x = 0. 1 ) 
and Moot* = 0.2) for several system sizes are presented in 
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FIG. 5: (a) System size dependence of the orbital correlation function 
Mqo (x) at x = 0. 1 , and (b) that at x = 0.2. 



Fig. [5] In Fig. [3a) foix = 0.1, Mqq(x) for several sizes cross 
around T/J = 0.25 below which Moo(^) increases with L. 
On the other hand, In Fig. [2tb) for x = 0.2, Moo(x) mono- 
tonically decreases with L in all temperature range. This dif- 
ference above and below x = 0.15 is also seen in the results 
of the correlation length. In Fig. [6] a correlation length at 
x = 0.15 and x = 0.2 are compared. In x = 0.15, %(x) for 
different sizes cross around T /J = 0.25. As shown in the in- 
set of Fig. |6ja), the scaling analyses works well. From this 
analyses for £ (jc), the OO temperature in x = 0.15 is obtained 
as 7bo(jc = 0.15)// = 0.248 ±0.003. On the other hand, in 
x = 0.2 [see Fig. |3fb)], £, (x) for different sizes do not seem 
to cross with each other at a certain temperature, and the scal- 
ing analyses does not work. From the above numerical re- 
sults, it is thought that the long-range OO disappears around 
0.15 <x< 0.2. 

The impurity concentration x dependence of Tqo(x) ob- 
tained by the MC and CE methods are presented in Fig. [7] 
Two kinds of the CE methods, where the PS operators are 
treated as classical vectors and quantum operators, are carried 
out. These are termed the classical and quantum CE meth- 
ods, respectively. In both cases, we adopt the two-size cluster. 
As a comparison, the Neel temperatures in the 3D XY model 
obtained by the classical MC method, and those in the 3D 
Heisenberg model by the classical CE method are also plot- 
ted in the same figure. It is shown that decrease of Too(x) by 
the MC method is much steeper than that of 7^ (x) in the XY 
and Heisenberg models. As shown in the size dependences 
of Moo(x) and §(x) at x = 0.2, it is thought that the long 
range OO is not realized at this impurity concentration. A 




FIG. 6: (a) System size dependence of the correlation length £ (x) /L 
at x=0. 15, and (b) that at x=0.2. The inset of (a) is the scaling plot 
for £, (x) at x = 0. 1 . The OO temperature and the critical exponent at 
x = 0.15 are obtained to be T 00 (x) = 0.248 ±0.003 and v = 0.755 ± 
0.085, respectively. 
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FIG. 7: Impurity concentration x dependence of the OO temperature 
Tqq(x). Filled circles are obtained by the MC method. Results by 
the quantum and classical CE method are shown by broken lines. 
For comparison, x dependence of the Neel temperature T^(x) in the 
3D XY model obtained by the MC method and that in 3D Heisen- 
berg model by the classical CE one are presented by filled triangles 
and dotted line, respectively. Thick arrow indicates the percolation 
threshold in a 3D simple cubic lattice. 
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FIG. 8: (a) A snapshot in the MC simulation for the PS configuration 
at x=0. 1, and (b) that at x = 0.3. Filled circles indicate impurities. 



rapid decrease of Too (x) in comparison with the spin order- 
ing temperatures is also obtained by the CE method. The OO 
temperature monotonically decreases with x, and disappears 
around x = 0.4 in the quantum CE calculation, and around 0.5 
in the classical CE one. The critical impurity concentrations 
obtained by the MC and CE methods are much smaller than 
the percolation threshold x p = 0.69 in the 3D simple-cubic 
lattice. 



Let us explain the physical picture of the orbital dilution. 
Snapshots of the PS configuration in the MC simulation are 
shown in Figs. [8ja) and (b) for x = 0.1 and 0.3, respectively. 
The staggered-type OO with the orbital angle (6^, 6b)=(Q,ti) 
is seen in the background of Fig.|HJ a )- At the neighboring sites 
of the impurities indicated by the open circles, PS vectors tilt 
from the angle of (0, n). This deviation of the PS angles is 
not only due to the thermal fluctuation. Focus on the NN sites 
along the x direction of an impurity which occupies the down 
PS sublattice. In almost all these sites, PS angles are changed 
from to a positive angle 80. This kind of tilting from (0, n) 
becomes remarkable at x = 0.3. Then, we explain the micro- 
scopic mechanism of this PS tilting due to dilution [see Fig.[9]|. 
Focus on a PS at a certain site termed i. The interaction act- 
ing on this site is considered by the MF approximation where 
we assume the staggered-type OO with the PS angle (0, n) ex- 
cept for the site i and an impurity site. The Hamiltonian which 




■x V 



FIG. 9: (a) A schematic PS configuration without impurity, and (b) 
that with an impurity. A filled circle represents an impurity. 



concerns the interaction acting on this site is given as 



= - I hi Ti 

/=(-w) 



■ £i-e, Tf_ 



«/ H-e, / "i 



(24) 



where e/ is a unit vector along / in the simple cubic lat- 
tice, and h; = (h*,^) are the MF. In the case of no 
dilution [Fig. |9|a)], the mean-fields are given by h x = 
/(-V^/2,-1/2), h y = 70/3/2,-1/2) h- = 7(0,-2), and 
the Hamiltonian in Eq. ( f24b is reduced to 



v/ 

'I 7 



(0 



3/7r. 



(25) 



This implies that the stable PS configuration at the site i is 
0j = K. Then, introduce an impurity at site ; — e x and consider 
the PS at site i [Fig. [9jb)] . The ^-component of the MF in the 
case without impurity is changed into h x = J(— \/3/4, — 1/4), 
and others are not. The effective interaction in Eq. (l24l is 
given as 



(0 



= J -[ 



(26) 



implying that the stable orbital angle at site i is 0,- ~ % — 0.15. 
This PS tilting due to dilution is attributed to the fact that 
the orbital interaction explicitly depends on the bond direc- 
tion and is the essence of the diluted orbital systems. This 
is highly in contrast to the dilute spin system where dilution 
does not cause specific spin tilting around the impurity site but 
simply increases thermal spin fluctuation since number of the 
interacting bond is reduced. 



V. DILUTION IN THE SPIN-ORBITAL MODEL 

In this section, we examine the dilution effect in the spin- 
orbital coupled model described by Mst m Eq. ©. First, 
we briefly introduce the MF calculation for the spin and or- 
bital structures at x = 0. The two sublattice structures for 
both the spin and orbital ordered states are considered, and 
the PS angles in sublattices A and B are assumed to be 
(6a, Ob) = (6,— 9). We obtain the ferromagnetic spin order in 
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FIG. 10: (a) Total energy E, and (b) the A-type AFM correlation 
function Ma- af ( x = 0) calculated for several values of the minimum 
energy E m \ n in the WL method. System size is chosen to be L = 10. 



the case of J1/J2 > 3, and the A-type AFM one in J1/J2 < 3. 
In the A-type AFM state, the orbital PS angle is uniquely de- 
termined as 6 = cos {2/2/ (5/i — J2 + 6g)}. By taking the 
MF results into account, for the following MC calculations, 
we choose the parameter set as (J\j Jiigj Ji) = (2.9,5). In 
these values, the OO appears at much higher temperature than 
the Neel one, and the A-type AFM is realized near the phase 
boundary between FM and A-type AFM. These are suitable 
to demonstrate the magnetic structure change due to dilution. 
The MC simulation results in the realistic parameter set for 
LaMn03 will be introduced in the Sect. [Vl] 

In the MC simulation, we utilize the WL method in L x L x 
L site cluster (L = 6 — 10) with the periodic -boundary condi- 
tion. The spin operator Sj in the Hamiltonian is treated as a 3D 
classical vector with an amplitude of 1/2. In the simulation, 
2x 10 7 MC steps are spent for measurement after calculating 
the histogram for the density of states. Physical quantities are 
averaged over 10MC samples at each parameter set. We no- 
tice again the lowest energy edge E m i n in the density of states 
which is introduced in Sec. III. In Fig. [10] we show the E m i n 
dependence of the total energy, E, and A-type AFM correla- 
tion function, Ma-afM defined by 



M 



A-AF 



N{\ -x) 




1/2 



-1)''*S, 



(27) 



where // for / = (x,y,z) represents the I component of the co- 
ordinate at site i. The results in Fig. fTUl a) imply that the tem- 
perature below which the total energy is flat is determined 




FIG. 1 1 : Temperature dependence of the orbital correlation function 
Moo (x = 0) , the PS angle function M ana (x = 0) , the A-type AFM one 
A/a-Af( x = 0), and the FM one Mp(x = 0). System size is chosen to 
be L = 8. 



by an adopted value of iSmin- This temperature is denoted 
as Tmin from now on. As shown in Fig. [TOfb), in the case 
of E min = -9.75 (-9.84), an obtained M A _ AF (x = 0) below 
Jmin is about 55% (75%) of its maximum value of 1/2. That 
i s ' 7mi n at Emj n = —9.84 is lower than the Neel temperature of 
the A-type AFM. Although a saturated value of M A _ AF (x = 0) 
is less than 0.5, this result is enough to examine the ordering 
temperature. We chose E m i n = —9.84 in the following MC 
simulation and focus on change of the magnetic ordering tem- 
perature due to dilution. 

First, we show the results without impurities. In Fig. [TT1 
calculated Moo (x = 0), M A _ A p(x = 0), and the ferromagnetic 
correlation function defined by 



M F (x) 




L/2 



(28) 



are presented. The staggered-type orbital correlation func- 
tionMooO* = 0) abruptly increases around T jJ% = 2.5 which 
corresponds to the OO temperature Tqq{x = 0). This value 
is consistent with the previous results obtained in the model 
Hamiltonian the effective orbital interaction in the 

present Hamiltonian ,3^st with paramagnetic state is / or t, = 
g + 3/i/4 — J2/4 where S; • Sj in J%st is replaced by zero. 
The obtained Too{x = 0) = 2.572 corresponds to 0.3/ or b in 
the present parameter set. This value is consistent with 
r 00 (> = 0) = 0.3447 obtained in Sect. IV [see Fig. Ha)]. 
In Fig. [TTJ the angle correlation function M ang (x = 0) starts 
to increase at Tqq(x = 0). With decreasing temperature, at 
around T/J2 = 0.5[= T^(x = 0)], the second transition oc- 
curs. The orbital correlation function MooG* = 0) decreases 
abruptly, and Ma-afG* = 0) grows up. The ferromagnetic 
correlation function Mp(x = 0) shows a small hump structure 
around 7n(x = 0). That is, T^(x = 0) is the Neel tempera- 
ture of A-type AFM. The PS angle correlation M ang (jr = 0) 
decreases and almost becomes zero below 7n(x = 0). This 
result indicates that, due to the magnetic transition, the PS 
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FIG. 12: (a) System size dependence of the orbital correlation func- 
tion Moo ( x = 0), and (b) that of the A-type AFM correlation function 

M A -af(* = 0). 



FIG. 13: (a) Impurity concentration x dependence of the A-type 
AFM correlation function Ma-afM, and (b) that of the FM cor- 
relation function Mp(x). System size is chosen to be L = 8. 



angle is changed into (6a, 0b) ~ (tt/2, — 7T/2) which is con- 
sistent with the MF results. In Figs. [12] size dependences of 
Ma-af{x = 0), MooC* = 0) are presented. With increasing 
L, changes of Mqq(x = 0) and Ma-af(x = 0) at 7n(jc = 0) 
become steep, although a saturated values of Ma-af( x = 0) 
is still less than 0.5 due to a finite value of |.E m j n — Eqs\ as 
mentioned above. 

Impurity concentration x dependences of Ma-af(x) and 
Mp(x) are presented in Fig. [T3] With increasing x from the 
x = case, Ma-af( x ) decreases gradually and almost dis- 
appears around x = 0.09. On the other hand, Mp(x), which 
shows a small hump structure around T/J2 = 0.4 at x = 0, in- 
creases with x, and takes about 0.35 in the case of x > 0.09. 
That is to say, the magnetic structure is changed from A- AFM 
into FM by dilution. At x = 0.06, both Ma-af (x) and M F (x) 
coexist down to the lowest temperature in the present simula- 
tion. This is supposed to be a cant-type magnetic order or a 
magnetic phase separation of the FM and A-type AFM phases. 

To clarify the mechanism of the magnetic structure change 
due to dilution, the effective magnetic interaction and the PS 
configuration are examined. Here, the AFM stacking in the 
A-type AFM structure is chosen to be parallel to the z axis. 
The effective magnetic interaction j\ is defined such that the 
Hamiltonian Mst in Eq. (0 is rewritten as = Jj^i ' 
Sj. The explicit form of the effective interaction is given as 



jf = 2(j l +j 2 ) 



2J 2 lT? + Tf 



3 1 



(29) 



where we consider a NN pair of sites ; and j(= i + e,) along 
the z direction, since we are interested in the magnetic struc- 
ture along z- A contour map of the effective interaction Jf, 



and a snapshot of the PS configurations in the same xy plane 
are presented in Fig. \l4\ a) and (b), respectively. Signs of Jf 
in almost all region are positive (antiferromagnetic), reflecting 
the A- AFM structure. At the neighboring sites of the impurity 
along the y direction, Jfs are negative (ferromagnetic). Away 
from the impurity, PS are ordered as ±T X in the staggered- 
type OO. Near the impurity, PS tilt from ±T X and finite com- 
ponents of T~ appears. This tilting of PS is seen in the results 
of J%t as explained in Sec. [IV] Based on these numerical 
simulation, we explain mechanism of the magnetic structure 
change due to dilution. Start from the staggered-type orbital 
ordered state of (T x , T z ) = (±1/2,0). Introduce one impurity 
at a site z'o which belongs to the T x =1/2 sublattice, and focus 
on the PS configuration and the effective exchange interaction 
at sites i'o + e m and i'o + <?,„ + e z for m = (x,y) (see Fig. [T5l ). 
As explained in Sect. IIVI orbital dilution induces the PS tilt- 
ing so as to gain the energies of the bonds where an impurity 
does not occupy. Thus, PS at site z'o + me tilts from 9 = 3n/2 
to 37r/2 + 89(— 80) for m =x(y) with a positive angle 89. 
Since the orbital interaction is the staggered-type, the bilin- 
ear term T, z , . 7? , - , - in Eq. ( |29| > are negative for both the 
m = x and y cases. As for the linear term in Eq. ( |29l , there is 
a relation (7^ + 7*^) = -(7*^ + 7*^^). That is, 
contribution of this linear term to the spin alignment along z, 
which is determined by a sum of J 7 , , and Jf , - , is canceled 

J 'O+^.v 'Q+ e y 

out. Therefore, when the first term in Eq. d29l overcomes the 
positive constant 3J 2 /2 —J\/2, Jf becomes negative and the 
ferromagnetic alignment along the z direction is stable around 
the impurity sites. 
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FIG. 14: (a) Contour map of defined in Eq. j29\ , and (b) a snap- 
shot of the PS configuration around an impurity in the xy plane ob- 
tained in the MC method. A filled circle represents an impurity. Tem- 
perature is chosen to be T/J2 = 0.3. 



VI. SUMMARY AND DISCUSSION 

In this section, we discuss implications of the present nu- 
merical calculations on the recent experimental results in the 
transition-metal compounds. First we have remarks on the re- 
lation between the calculated results of J%t shown in Sect.lIVI 
and the experiments in KCui_ v Zn v F3. 10 As shown in Sect llVl 
Too i x ) rapidly decreases with increasing x in comparison with 
dilute magnets (see Fig. |7). Although the critical concentra- 
tion (x ~ 0.2 — 0.5), where the OO disappears, depends on 
the calculation methods, that is, MC and CE, these values are 
far below the percolation threshold (x p = 0.69). This result is 
consistent qualitatively with the Zn concentration dependence 
of the OO temperature in KCui_ x Zn x F3 where OO vanishes 
around x = 0.45. One of the discrepancies between the the- 
ory and the experiments are seen in their quantitative values 
of the critical impurity concentration where OO disappears. 
Some of the reasons of this discrepancy may be attributed to 
the anharmonic JT coupling and the long-range PS interac- 
tions due to the spring constants beyond the NN ions and so 
on, both of which are not taken into account in the present cal- 
culation. The former subject, i.e. the anharmonic JT coupling, 



P 
z 




y T 



FIG. 15: A schematic PS configuration around an impurity at site ;'o- 
A filled circle represents an impurity. 



induces the anisotropy in a bottom of the adiabatic potential 
of the Q x — Q z plane, and prevents the PS tilting around impu- 
rity sites. This effect on the reduction of Too(x) was studied 
briefly in Ref. [HI It was shown that, in the realistic parame- 
ter values, the reduction of Tqo(x) becomes moderate by the 
anharmonic coupling, but it is still steeper than that in dilute 
magnets. Another factor which may explains the discrepancy 
between the theory and the experiments is the quantum aspect 
for the orbital degree of freedom. In the results obtained by 
the quantum CE method as shown in Fig. [7] the critical x for 
Too( x ) = is larger than the results by other two classical 
calculations for the orbital model and is close to the experi- 
mental value of x = 0.45. This may be due to the fact that 
quantum fluctuation in low temperatures weakens the low di- 
mensional character in the OO state and prevents a collapse of 
OO against dilution. This kind of quantum effects in the dilute 
orbital system was examined by the present authors in the two 
dimensional quantum orbital model 3 It was shown that the 
reduction of Too(x) due to dilution is weaker than that in the 
classical orbital model. 

We briefly mention the orbital PS tilting due to dilution. 
Similar phenomena are known as a quadrupolar glass state in 
molecular crystals where different kind interactions between 
molecules with quadruple moment coexists^ A kind of glass 
state in terms of the quadrupole moment appears with increas- 
ing randomness for the interactions. We suggest a possibility 
that the present observed PS tilting accompanied with the lat- 
tice distortion of ligand ions is able to be detected experimen- 
tally. One of the most adequate experimental techniques are 
the pair-distribution function method by the neutron diffrac- 
tion experiments, and X-ray absorption fine structure (XAFS) 
where the incident x-ray energy is tuned at the absorption edge 
of the impurity ions. This observation may work as a check 
for the present scenario in the dilute orbital system. 

Next we discuss implications of the calculated 
results in Sect. [V] to the experimental results in 
LaMni_^Ga r Q 31 13 i 14 i 15 i 16 - 17 - 18 By analyzing the spin- 
orbital coupled Hamiltonian J%st, we find that the magnetic 
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FIG. 16: Impurity concentration dependence of the A-type AFM 
transition temperature, 7n(x), and the FM one, T c (x), calculated in 
the realistic parameter values for LaMni_ v Ga v 03. The parameters 
and the system size are chosen to be (Ji/Jiig/Jl) — (2-5,5) and 
L = 8, respectively. 



structure is changed from the A-AFM order into the FM 
one. This calculation qualitatively explains the experimental 
results in LaMni_ Y Ga r 03 from the macroscopic point of 
view. In Sect. [V] the parameter set is chosen to be close to 
the values for the A-type AFM/FM phase boundary, in order 
to demonstrate clearly the magnetic structure change due to 
the orbital dilution. Here we briefly introduce the numerical 
results obtained in the realistic parameter values. To evaluate 
the realistic values, we calculate the OO temperature, the 
Neel temperature by the MF approximation, and the spin 
wave stiffness by the spin wave approximation from Jtj$r, 
and compare the experimental results in LaMn03. Then, 
we set up the parameters as (/i/J^g/jy = (2.5,5). The 
x dependences of the magnetic transition temperatures are 
presented in Fig. [16] With increasing x from x = 0, T^(x) 
of the A-type AFM order gradually decreases, and around 
x = 0.2, the A-type AFM is changed into the FM order 
which remains at least to x = 0.4. In semi-quantitative 
sense, this result is consistent with the experimental magnetic 
phase diagram in LaMni_ Y Ga Y 03. However, one of the 
discrepancies is that the canted phase survives up to x = 0.4 



in LaMni_ x Ga v 03. This difference between the theory and 
the experiments is supposed to be due to the t2 g spins in Mn 
sites and the antiferromagnetic superexchange interaction 
between them which are not included explicitly in the present 
calculation. This interaction stabilizes the A-type AFM phase 
in comparison with the FM one, and maintains the canted 
phase up to a higher x region. 

In summary, we present a microscopic theory of dilution ef- 
fects in the e g orbital degenerate system. We analyze the dilu- 
tion effects in the e g -orbital Hamiltonian without spin degree 
of freedom, J#r, and the spin and e g orbital coupled Hamilto- 
nian, J#st- The classical MC simulation and the CE method 
are utilized. It is shown that the OO temperature decreases 
rapidly with increasing dilution. From the system size depen- 
dence of the orbital correlation function in the MC method, 
the OO is not realized at the impurity concentration x = 0.2. 
Tilting of orbital PS around impurity is responsible for this 
characteristic reduction of Too(x). This is consequence of 
the bond dependent interaction between the inter-site orbitals. 
In the analyses of the spin-orbital coupled model, the mag- 
netic structure is changed from the A-type AFM structure into 
the FM one by dilution. This is explained by changing of 
the magnetic interaction due to the orbital PS tilting around 
the impurity. The present results explain microscopically the 
novel dilution effects in KCui_ jr Zn JC F3 and LaMni_ Y Ga Y 03, 
and provide a unified picture for the dilution effect in the or- 
bital ordered system. 
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